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ABSTRACT 


This work presents a method for obtaining optimal topology for stiffened plate 
structure. The method is based on a combination of existing finite element methods 
and evolutionary structural optimisation methods. Once the initial finite element 
problem is defined, the method proceeds by finding recovered smoothened stress field 
througout the structure and the elements are subdivided into subelements inorder 
to improve the accuracy of identification of the low stressed material. Low stressed 
regions are identified in the stiffened layer and remeshing is done seperately in each 
regions. The elements in the low stressed regions are removed by decreasing the 
material properties to near zero and this evolved structure is solved. The process is 
repeated by increasing the lower stress limit for plotting low stressed regions in each 
iteration, until the stopping criterion is reached. The results produced are presented 
and its advantages are highlighted. 
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Chapter 1 


INTRODUCTION 


1.1 GENERAL 

Rapid progress in computer hardware has opened new vistas for structural analysis 
and design. Large-scale problems with high requirements of in-core and out-of- 
core memory and processing speed have become tractable even on a PC. One such 
problem is that of structural optimization. This has been a topic of interest for many 
years, but has seen a recent resurgence with the rapid progress made by computer 
hardware. Typically, structural optimization was restricted to truss or frame type 
structures and laminated composites (lay-up, cut-out shape optimization). 

Topology optimization is a recent research topic, with a goal to start from a 
“slab” of material and “carve” out the final optimal design from it. This can lead 
to opening of new holes (i.e changing the topology) and “shaping” of existing holes 
in a structure, to lead to a fully stressed structure [4]. For example, let us consider 
the beam given in Fig. 1.1 (a) and Fig. 1.1(b). Starting from a solid beam, with the 
loading shown, material was removed iteratively to reach the final hollow beam 
configuration shown. Here it was required that the beam does not buckle under 
the applied (constant) axial load, and the Mises stress at any point does not exceed 
eighty percent of the yield stress. From the figure we note that the amount of 
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Figure 1.1: Topology optimization, of a solid beam(all dimensions are in meters) 


2 



material removal (and profile of hollow regions) depends on the applied loading, the 
end conditions and the applied constraints. In this case, the constraint on buckling 
dominated. 

In engineering components, plate or shell type structures are commonly used, 
e.g. wing panels, ribs, fuselage, cockpit of aircraft; bonnet and body of a car. Most 
of these structures are made of light material and “stiflfened” using longitudinal and 
transverse members (spars, ribs,longerons,bulk heads). For a typical aircraft, a sig- 
nificant amount of the total structural weight is due to the stiflfening components. 
Motivated by the simple example considered above, one can pose the following ques- 
tions: 

Is it possible to design the “optimal” stiffener configuration for a plate or shell-type 
structure? 

In principle, topology optimization should lead to the best possible (minimum 
weight) configuration. Available topology optimisation procedure can be classified 
into the following groups: 

a) Homogenisation based [14] 

b) Soft kill/Hard kill [6] 

c) Adaptivity based [3] 


1.1.1 Homogenization based topology optimization 

The domain is broken into a fine mesh of periodic cells(or subdomains) with the help 
of a very fine finite element mesh. Based on the current analysis of stress levels, 
due to an applied load, material is removed from a cell (element) whose Mises stress 
is below a lower threshold level. This may lead to a “porous” structure with filled 
and empty cells spread out. Getting a finite element solution of sufch fine mesh 
is impossible, and results may be useless due to spurious corners created in the 
domain(at vertices of hollow cells). Thus, the porous domain is homogenised and 
an average(smoothened) stress field is obtained, which is needed to propagate the 
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material removal(see [14]- [15] for details). 


1.1.2 Soft kill/Hard kill 

A fine mesh is made over the domain and material is removed(Hard kill) or replaced 
with a soft material(Soft kill) in elements with low average Mises stress(see [6]- [7] 
for details) . The new material domain is again solved for till convergence. 

1.1.3 Adaptivity based methods 

The goal here is to accurately capture the low stress regions by refining the mesh 
adaptively in regions of low stress. This leads to very small sized elements at the 
boundaries of the low stress regions (i.e. stress below a given lower threshold). Mate- 
rial is then removed from the region, and the process is continued till convergence(see 
[3]- [5] for details). 

The methods mentioned above suffer from the following drawbacks; 

1. Spurious corners are cieated in the domains, where material is removed from an 
element. This leads to stress singularities, which causes the process to go back and 
forth without convergence. This is because the regions from which material gets 
removed demands material to be brought back at the next step. 

2. Since average stress are taken, stress profiles are not very smooth. 

3. The homogenization procedure can give stresses that are very different from the 
initial stress state(see [2]), and hence wrong design. 

4. A fine mesh has to be resorted to in all cases, leading to heavy computational 
effort. 

The flaws of the existing methods can be circumvented, if the following goals can 
be met: 

1. Use of coarse mesh for the computation to reduce cost. 

2. Get rid of spurious stress concentration induced due to the material removal 
strategy employed. 
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The current study is an attempt to systematically develop a procedure for topology 
optimization, which can meet the goals given above. The procedure given in this 
study utilizes the following tools: 

1. Use a coarse mesh, with a good stress recovery algorithms, which should (in prin- 
ciple) give a better and smoothened stress field, as compared to the one obtained 
directly from the finite element solutions. 

2. Capability to plot smooth contours of stress levels corresponding to a given 
threshold stress level. 

3. Capability to remesh the new domains(with holes), with desired level of refine- 
ment at the interior boundaries. 

Obviously, accurate representation of the pointwise stress is very important to en- 
sure reliability of the design obtained. The recovered smoothened stress field can 
also be employed to get an estimate of the error in the approximation, which can 
be used for adaptive refinement of the finite element mesh obtained by remeshing. 


1.2 ORGANISATION OF THE THESIS 

The study in the thesis is organised as given below: 

1. The first chapter deals with a discussion on the state-of-the-art and the goals 
of this study. 

2. In chapter-2, a higher order plate theory(HSDT) is given. This plate theory 
will be used to study unsymmetric laminated plates. 

3. In chapter-3, the finite element formulation corresponding to the HSDT model 
will be developed. 

4. Chapter-4 deals with the a posteriori recovery of a smoothened stress field, from 
that obtained using the finite element method. 

5. In chapter-5, the problem of an aluminium plate with an attached (underly 
ing) steel plate is considered. The goal being to remove material from the steel 
plate, inorder to obtain a light weight design, appropriately stressed and safe 
plate design. The subsections in this chapter describes how contours are plotted 
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for the low stressed regions from the recovered smoothened stress field. These 
contours are drawn using a grid of “subelements” embeded in an element. The 
contours are polygonal in nature (obtained by connecting adjacent stress field 
points with a straight line). A spline based curve is fitted to the polygonal 
contours, to get smooth enclosed regions of low stress. 

7. In the seventh chapter numerical problems are solved and results are shown. 

8. Finally, chapter eight gives the scope of this work followed by discussion on the 
future scope. 
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Chapter 2 


HIGHER ORDER PLATE 
THEORY 


The relationship between the strains at any point within a laminate and the 
corresponding deformations are the functions of the assumed displacement fields. 
The displacement field and strain-displacement relationships are given below: 


2.1 DEFINITION OF DISPLACEMENT 
FIELDS 

The membrane-flexure coupling phenomenon exhibited by an unsymmetric lami- 
nate, necessiates the use of a displacement field containing both. The displacement 
field Yix,y,z) — [u{x,y,z), v{x,y,z), w{x,y,z) ] is derived from the expanded 
Taylor’s series in terms of thickness coordinate {z) is considered. This is given as, 
[ 1 ] 

u{x, y, z) = Uq{x, y) -h z9^{x, y) -I- z‘^(t>x{x, y) + z^'ijj^{x, y)\ 

v{x, y, z) = vq{x, y) -f zey{x, y) -I- z^(l>y{x, y) + z^-ijjyix, y)\ (2.1) 

w{x,y,z) = wo{x,y) 
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In the expansion (2.1), it is assumed that transverse normal strain is zero. 

The linear strain-displacement relationships using small deformation theory can be 
written as follows: 

^xx ^0)1 “b Z Oxix ~b ^ ”b ^ V^aciD 

^yy = V0,y + 2 : ey,y + (f>y,y + Ipy.y, 

Ixy = «0,y + -2 Qx,y + + 2® V’x.j, 2) 

-bUO,,, + 2 Oy,^ -h (t)y,^ -f- 

'Yyz “ Oy -|- 2 .Z (f)y 3 Z Ipy ”1“ U^0)y) 

Tx2 ” ^x "b 2 Z (j)x “b 3 Z ll)x "b U^Oja: 
where comma (,) denotes the partial derivatives. 



Figure 2.1: Representation of arbitrary three dimensional domain 

The condition that the transverse shear stresses vanish on the plate’s top and bottom 
faces (see Fig. 2.1) is equivalent to the requirement that the corresponding strains 
be zero on these surfaces, i.e. 


7yz{x, y, ±|) = Jxz{x, y, ±|) = 0 


On introduction of the conditions as given above in the expressions for transverse 
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shear strains, the following relations are obtained. 


= = 0 and 


(2.3) 


'^y — 3/i 2 (^y "h '^^Oiy) ) '4’x 3/i2 (^x "b ‘^Otx) 

The displacement field of Eq. (2.1) is modified by setting 4>x and to be zero 
according to conditions of Eq. (2.3). The resulting displacement field is written 
below; 


u{x, y, z) = uo(x, y)+z y) + ipx{x, y)\ 

v{x, y, z) = vq{x, y) + z 6y{x, y) + z^ ipyix, y); (2-4) 

w{x,y,z) = woix,y) 


In Eq. (2.4) u, v and w are the displacements along x, y and z directions respectively. 
«o, Vo and wq are the mid-plane displacements while 9x, By are rotations about y and x 
axes, respectively. Whereas, V’x and V’y are higher order terms in the Taylor’s series 
expansion and are also defined at mid-plane. Thus, the generalised displacement 
vector {5} of the mid-surface contains seven degrees of freedom (DOF) and is given 
by; 


■{ 5 } — {^^ 0 , VQj Wo, dx, 9y, 'Ipx, 

Corresponding strain-displacement relationship becomes; 

fix ” '^051 “b ^ ^xjx "b ^ V’xjxi 
^yy ~ ^O’y ^ ^y>y 4" ^ ^y’y> 

Ifxy ~ ^0)y "b Z Bxiy *b Z 'ij^xjy 

-^Vo,x + Z 9y,^ -b z® 'lpy,j.-, 

lyz =9y+Sz‘^^y+ W0,y; 

7xz -9x + ^z^rpx + Wo,x 
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2.2 LAMINATE CONSTITUTIVE 


EQUATIONS 


A unidirectional fibre reinforced lamina is treated as an orthotropic material 
whose material symmetry planes are parallel and transverse to the fiber direction. 



Figure 2.2: Co-ordinate axes in lamina 

Following fig. (2.2) the material co-ordinates axes L and T are defined parallel 
and perpendicular to the fibre direction respectively, while global co-ordinates are 
X and y. The angle between global co-ordinate axis x and fibre direction, material 
co-ordinate axis L, is a and is known as the orientation angle. As a sign convention 
anticlockwise angle is taken as positive. 

Generalized Hooke’s law 

In the formulation of lamina constitutive equation following two assumptions are 
made. 

1. The lamina is a continuum. 

2. It behaves as a linearly elastic material. 

Further, at the micro-level the following assumptions are made about the material: 

1. Perfect bonding between fibres and matrix exists. 

2. Fibres are parallel and uniformly distributed throughout. 
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3. The matrix is free of voids or micro cracks and initially in a stress free state. 

4. Both fibres and matrix are isotropic and obey Hook’s law. 

Stress-strain relations for the lamina in the material coordinate axis, whose fibers 
are oriented at an angle a with reference to the x axis is given as ( [?], [?]): 


I (Jj — [Qtj]i I Cj (2-6) 

where {ctj} is the vector of stress components, [Qtj] is the stiffness matrix, and {e^} 
are the engineering strain components, for the lamina. 

The stresses and strains in the x, y and z directions are obtained by transformation 
relations given in equation (2.6). The transformed stress strain relations for the Z*^ 
lamina are given as: 


✓ N 

X 


Qii 

Qi2 

Qi6 

0 

0 



Gy 


Q21 

Q22 

Q26 

0 

0 


Cy 

^ '^xy 

^ = 

Qi6 

Q26 

Qee 

0 

0 

< 

Ixy ^ 



0 

0 

0 

Q44 

<345 


lyz 

^xz 

1 

0 

0 

0 

Q45 

<355 

1 

Ixz 

V. y 


(2.7) 


Using the above lamina constitutive equations and integrating the stresses over the 
laminate thickness, the stress resultants in terms of inplane forces, moments and 
shear forces per unit length, for a laminate with NL laminae are. 


' {N} ' 


' {^ 0 } ’ 

{M} 


p„] [DJ 0 


{«} 

{M*} 

^ = 

[OJ’’ [D,l 0 

i 

{k*} > 

{Q} 


0 0 [D,] 







( 2 . 8 ) 


where {eo} are the mid-plane strains, {«} are the mid-surface curvatures and {«*}, 
{i/}, {i^*} are higher order terms and 


rzi-i 

{N}’' = (W.,iV„iV^)=E/ (2-9) 

1=1 
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( 2 . 10 ) 


rzi-i 

~ (-Wa;, ^ ) *^X3/) ^ dz^ 

i=l 

{M*}^ = (M:,M;,M4) = / (cr,,a„a,,)' d2); (2.11) 

/=i 

pZl — i 

{Q}^ = (Qx, Qy) = E / (^-> (2.12) 

i=l 

rzi-1 

{Q*}T = (Q*, Q*y)=Yl [ {r^z. ^yzY dz. (2.13) 

1=1 

Here, z^ are the 2 : coordinates corresponding to the laminae interfaces as shown in 
Fig 2.3. 

The other terms and following rigidity matrices are given in Appendix A. 


Middle Plane 
X 


Figure 2.3: Geometry of multilayer laminate 

Thus, with the assumed displacement model, the various rigidity matrices derived 
are: 

[Dm] = Membrane; [Dj = Membrane-flexure coupling. 

[Db] = Flexure; [Dg] = Shear. 

The set of rigidity matrices [Dm], [Dc], [Db] and [Dg] are used in forming overall 
rigidity matrix [Dr] for the laminate. 
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Chapter 3 


FINITE ELEMENT 
FORMULATION 


Triangular elements are used for finite element discretisation and hierarchic shape 
functions of order p (p < 4) are used for finite element approximation. The mesh 
generation is done using automatic mesh generator and a mesh generated over the 
plate domain is shown in fig. (3.1). 



Figure 3.1: Mesh generated over a square domain 




3.1 DEFINITIONS 


Let {V} be the displacement vector defined as: 

( ^ 
u 

{V} = \ V [ (3.1) 

w 

Let the stress and strain vectors corresponding to {F} be {cr} and {e} which can 
be defined as: 



Prom the generalized Hooke’s law which relates stress components to the respective 
strain components in global coordinate system 

M = [Q] {€} (3.3) 

where material stiflfness matrix [Q] for orthotropic material is as given in Eq. (2.7) 
for each lamina. 

In any elastic body, linear strain-displacement relationships using small deformation 
theory is given in Eq. (2.5). 

The above relations in the matrix form can be expressed as: 

{£} = [D] {V} (3.4) 


where [D] is a differential operator in terms of global coordinates such that 


— 0 
dx ^ 

^ dy 

A. A 

dy dx 
0 Tx 

.f. 0 



(3.5) 
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The 

in- 


+. be written in terms 

components of displacement can 

whichi can be written as. 
plane functions, which 


.tnr is written in matrix form as 

The displacement ve 

(V} = l^] 


of the seven unknown 


( 3 . 6 ) 


where 


1 0 0 z 6 ^ 

l$l=l 0 1 0 0 S 0 

0 0 1 0 0 0 0 



( 3 . 7 ) 


„„ etc are approximated 

. +Vip functions uo, ^o, - 

I, ,te finite element app—on ‘fie fit 

using nsfiape functions per elemen. 


( 3 . 8 ) 



W = t 

i=l 


( 3 . 9 ) 


generaltod displacement vector 

element. ^ independent coefficients 

• oUnne functions are usea, 

KeTnorfc Wfien Lagrangian 

correspond to pfiysiea' noiee. 
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Hence displacement vector can be written as 


{v'} = m [w] {d} (3.10) 

where {^2}^, and [N] is a matrix in terms of n inplane 

hierarchic shape functions given in Appendix B. 

Hence, the strain vector can be written as 

{e} = {[D] [$] [N]) {d} (3.11) 

and the stress vector 

{a} = [Q] ([^] [#] [iV]) {d} (3.12) 


3.2 FINITE ELEMENT FORMULATION USING 
ENERGY PRINCIPLE 


The total potential for the plate is given by 

n (.5 ) = Ejt* (5). (3.13) 

e=l 

where tt® is the total potential of the non-intersecting (but adjacent) sub-domians 
e which are part of the domain (N sub-domains are considered here). The total 
potential can be expressed in terms of internal strain energy and external work 
done kF(®), as follows; 

7r® ((5) = - W(®). (3.14) 

Strain energy of the laminate can expressed as follows: 

= 5 

Work done by the applied external transverse load is, 

(3 16) 

= Vw wo r wo f- 
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where /+ is the transverse load on the top face and f~ is the transverse load on 
the top face R~. 

The exact solution Uex to this problem is the minimization of the total potential II. 
This can be obtained as: 


n = 0, 

which is also the Virtual Work Formulation of the problem in terms of the com- 
ponents of 5. Prom this, the seven coupled equilibrium equations in terms of the 
components of 5 can be obtained. This leads to generalised finite element formula- 
tion: 

[K] {d} = {P} (3.17) 

where [K] is the stiffness matrix and {F} is load vector. 

In the next section, we are going to employ the virtual work formulation to derive 
the finite element formulation of this problem. 


3.3 COMPUTATION OF ELEMENT STIFFNESS 
MATRIX 

The stiffness matrix corresponding to assumed deformation state of an element 
can be defined by expressing the internal strain energy in terms of unknown nodal 
displacements. In the formulation of unsymmetric laminates the membrane, the 
flexure,membrane-fiexure coupling and shear strains contribute to strain energy. 
The set of Eq. (2.5) can be used along with Eq. (3.12), (3.15) to express these 
strains in terms of nodal displacements. 

The internal strain energy of an element (given by area A®) can be determined as: 

U‘ = I L. ( {dV ( [Sm]"' Pm] [Bml ) {<*} + ( PI"' [A] [A] ) {<«} 

+ {dY ( [B.f [A] PI ) {i} + {dY ( Pm]"' [A] [P] ) {<!} (3.1S) 

+ {dY ( P]"' [A] Pm] ) {d} ) dA 
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Here the first term on the right-hand side of Eq. (3.18) is the in-plane contribution, 
the second and third terms are out-of-plane bending and through-thickness shearing 
contributions, respectively. The last two terms are contributions from coupling 
between in-plane and out-of-plane actions. 

The strain energy expression can be written in a concise form as: 

= i {dy [K^] {d^} (3.19) 

where \K^] is the stiffness matrix and {d®} are the coefficients corresponding to 
the finite element solution in an element e. All the terms in equation (3.18) are 
evaluated individually and then summed to yield [AT®]. 

It is given by: 

[K’] = /^. ( ( P™! [B„l ) + ( [B»F [A] [Bd ) + ( [B,Y m [B.] ) 

+ ( [BJY [Dc] [Bj] ) + ( [BtF [A] 1B„1 )]dA 

The numerical integration is carried out to get element stiffness matrix. 

The matrices [Bm], [Bb] and [Be] and vector {d} are given in Appendix B. 


3.4 COMPUTATION OF ELEMENT LOAD 
VECTOR 

For the extension and/or bending problems of laminated plates, the applied 
external forces may consist of following cases of loading: 

1. Uniformly distributed load acting over the element in the 2 :-direction on top or 
bottom bounding planes of the plate. 

2. Sinusoidal distributed load acting over the element in the z-direction on top or 
bottom bounding planes of the plate. 

3. Traction load acting along the edges in the x or ^/-direction . 

The total external work done by these forces on an element can be expressed as 
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follows: 


or 


{df 1 ([«] iivj)^{r.} ds 

(3.21) 

= {df {F^} 

(3.22) 


where {i^®} is the element load vector, and the traction vector is given as: 


{Tr} = 


T, 

Ty 


> 


T. 


(3.23) 


Transverse loading, % can be sinusoidal as Tz = Pc sin ( 2 ^) sin(^^) or uniform 
pressure = Pq . Tx, Ty can be uniform inplane tractions along x and y, repectively. 


3.5 GEOMETRIC APPROXIMATION 


3.5.1 STRAIGHT EDGE ELEMENTS 


The geometry is expressed in terms of the shape functions as: 



Xt 

II 

[ J 


(3.24) 


where, 


Xi 


Vi 


> are the coordinates of the node of the element, and Ni are the 


linear shape functions of the element. 


LINEAR MAPPING; 


Linear mapping in general is used when all sides of mapped elements are straight 
lines(see Fig. 3.2). Consider the mapping of variables from x,y to rj such that 
{x, y G A} on substitution of iV) these can be written in the matrix form as 


X — Xi 

y-Vi 


X2 - Xi 
2/2 - 2/1 


X3 -Xi 

Vz-yi 



(3.25) 
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Figure 3.2: Linear mapping in two dimension 


Clearly, 


^=X2-Xi 

% = y2-yi 


= ^3 - 


yi 


On inverting the matrix it is easy to see that 



1 

1 

CO 

1 

1 

H 

1 

CO 

T 


X — Xi 

. ^ J 

' ” A 

_ -(j /2 - yi) 

{X2 - Xi) _ 


[ y-yi 


where, A = (x 2 - Xi) {yz - yi) - {xz - xi) (y 2 - yi) 


Hence, 


( 3 . 26 ) 


( 3 . 27 ) 


^ = y.3.zyi. 
dx A 

drj _ V2-yi 

dx A 


§i — xa-xi 
dy A 

dri _ X2—XI 
dy A 


( 3 . 28 ) 


Transforming dxdy into d^dr}: 

Let the differential area dxdy is formed through vectors dx and dy with magnitude 
dA and direction normal to the elemental area is k 


dA = dx dy = [dx x dy]-k 


( 3 . 29 ) 


dxdy = [i^^d(i + ^dr,j)x(^dit + ^dr,])]-k 
= 1 J| d^ dy 
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or the Jacobian matrix is given by 


dx 

d]i 



dx 

in 

- ^ 1 ? 

Q'n 


Substituting from equation (3.5) 


I J| = (X2 - xi) (ys - y\) - (X 3 - xi) (y2 - Vi) 


(3.31) 


(3.32) 


Hence /e / (x, y) dx dy may be written as 

f f (x, y) dx dy= [ f (^, ri) d^ dr} (3.33) 

Je Jem 


such that / (^, T}) d^ dr] = | J| / (x(^, 77 ), y(^, 77 )). 


Numerical Integration: 

To compute f (^, 77 ) d^ dr} Gaussian quadrature for two-dimentional integrals 
is used through which 



/ (C> V) didr}^'^ f { ^,, 77 . ) w, 1 J1 

i=l 


(3.34) 


where, n is the number of quadrature points, Wi is the weight at quadrature 
point (^i, 77 t) and the Jacobian \ J\ is the area of the element e. 

Order of numerical integration is greater than or equal to (p + l )/2 where, p is 
the order of the integrand. 

Shape Function Derivatives: 

Shape functions are defined in the master element i.e. Nt = Ni{^, 77 ). Its deriva- 
tives with respect to x and y can be given as 

N- = I dN, drt 

di dx ‘ dr] dx /Q QC\ 

AT. — Ml . ^ §Mx , §ii V • / 

dy ^ dr] dy 


Since shape functions are defined in terms of normalised coordinates ^ and 77 the 
terms ^ and ^ can be evaluated at any point in the domain. 
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Chapter 4 


RECOVERY OF SMOOTHENED 
STRESS FIELD 

4.1 RECOVERY OF STRAINS 

Recovered strain-field over an element are constructed using patch-recovery 
technique [13]. In this recovered strains are obtained using minimization of energy 
norm of error over the patch (see Fig. 4.1) around an element J. Thus, the problem 
can be posed as: 


find e* such that 

/ = — / f ( CpB — €*)• Q { gpB - ^ ) dz dA (4.1) 

^ JApatch •'2=-t 

is minimized. Here e* and gpp are recovered and finite element strain vectors and 
Q is material stiffness matrix, respectively. 

The finite element strains and material stiffeness properties are already available 
over the patch. The recovered strains are constructed using the same mathematical 
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Figure 4.1: Patch over element J 


form as that of the exact solution of the model i.e. 


c XX 

-*0 

^ XX 

+ 2 

4- 

-♦2 
^ XX 

f* 

^ yy 



^ yy 

+ 

4- 

-♦2 
^ yy 


^ ^ xy 

+ ^ 

+.z^ 

e*2 

xy 


= Ylz + 7* 

1 

yz 



— -L. 

1 XZ ' ^ 1 XZ 



(4.2) 


In-plane strains of the strain tensor are approximated by p*'* order polynomial 


and are given by 


,♦0 


(p+i)(p+2)/2 

E 

i=l 


„»0 

XX, t 


Qi 


(4.3) 


Similarly, the strain components can be 

approximated as above(eqn. 4.3). Whereas, the transverse shear strains of the 
strain tensor are approximated by (p -I- i)*'* order polynomial and are given by 


(p+2)(p+S)/2 

E (4-4) 

t=i 

Similarly, 7 *^^, 7 *°^, 7*L strain components can be approximated as above(eqn. 4.4). 
The minimization problem 4.1 can be restated in terms of the unknown coefficients. 
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p *0 p*l ^*0 ^*1 ^ p *2 ^*0 /- y *! 

^ xx^i^^ XX, ^ XX, 25^ 2 / 2 /, i’ yy^t^ 2/2/j^’ xy,2?^ xy,2» xy,v / y2,25 / yz,i^ / |/z,t5 / xz,2’ / xZjZ? 

7*L,i- leads to a matrix problem over each patch. The monomials (shape 

functions) Qz used in the representation 4.3, 4.4 are given by 


?i = 1 
Qa = 
qi = x^ 
gio = 
gi3 = x^y^ 
qi6 = x^ 
gi9 = x^y^ 

where 


11 

II 

q5 = xy 

q& = f 

q& = x^y 


qu = x^ 

gi2 = x^y 

qu = xy^ 

Cn 

II 

qn = x^y 

gi8 = x^y 

gao = xy^ 

gai = f 


\ 

^ 1 

^ ^ 

> 

/ 

/ 

1 


(4.5) 


(4.6) 


Here {Xc, Yc) represents centroid of the element J and P{x, y) is any point 
in the patch. Solution of the matrix problem over a patch gives the coefl&cients 
corresponding to the expansion over the patch. The components of the recovered 
strain field are then represented in terms of lagrangian shape function (of the same 
order as that used for the polynomial fitting). Thus, we obtain “nodal” values of 
the recovered strain components, for any element J (see Fig. 4.2). 


(p+i)(p+2)/2 

(4.7) 

^ xxi.^ji yj) 

t=i 

NN 


i=‘ 

(4.8) 


Here Nj represents lagrange shape function associated with ji** node in terms of 
normalised coordinates ^ and rj ,NN represents number of nodes, qi{xj,yj) repre- 
sents monomials at nodal coordinates,e*° 3 .(xj,t/j) represents the recovered strain 
components at nodal coordinates. Similarly all the recovered strain components 
are obtained for all the nodes in the element. This process is carried out through 
all the elements and recovered strain components are obtained for all the element 


nodes. 
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Figure 4.2: J** element in the patch is represented in terms of lagrange quadratic 
element 

4.1.1 Special Case: Adjoining dissimilar material domains 



Figure 4.3: Patch over element J with adjoining dissimilar material domains 

Consider a patch as shown in figure 4.3, Here the shaded elements have different 
material properties as compared with unshaded elements. In this case, while recov- 
ering strains over the element J , the elements with different material properties 
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when compared with element will not be considered. For the figure shown while 
recovering strains for element J in patch information the elements J9 ,J10 ,J11 J12 
and J1 will not be considered, since it is known that when material properties 
changes then strains are discontinous at the boundaries where material properties 
changes. 


4.2 SMOOTHENED STRESS FIELD 

The recovered smoothened strain components are obtained by averaging recov- 
ered strain components over each node to get single value per node.This ensures 
continuity of strains at element boundaries. For example, figure 4.4(A) shows a 
patch of quadratic Lagrange elements. It has 6 nodes per element. Here recovered 
strain components for each element nodes are found separately by constructing a 
patch over each element. Because of this, the nodal strains in element A and in ele- 
ment B are not continous at their interface, thus you can visualize this effect as if the 
two elements are getting seperated by a marginal distance (see figure 4.4(B) (i)). In 
order to get continuity at element interface(see figure 4.4(B) (ii)), the recovered strain 
components are averaged at the common nodes, which the two elements shares, i.e., 
the recovered strain components of the nodes 3 and 8 in elements A and B are taken 
and their values are replaced by its average value respectively. Look at the node 13, 
it is shared by elements A,B,C,D,E,F,G and H so the recovered strain components 
of the node 13 in all the elements are taken and replaced by its average value. This 
process is repeated for each and every node to ensure continuity of strains over the 
element boundaries. Then the recovered strains are obtained at desired laminate 
thickness using the relations given in equation 4.2 . 

The plate model that is considered guarantees continuity of inplane stresses 
and strains, thus the inplane stresses are obtained through laminate constitutive 
equations which are given in equation (2.7) , but in out of plane direction the 
plate model guarantees only strain continuity and not stress continuity. Hence the 
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Figure 4.4: Sharing of nodes over element patches 

transverse stress field are obtained from the equilibrium equations. 

To make engineering decisions one should guarantee that the finite element 
stresses that are obtained are sufiiciently close to exact solution of the mathematical 
model considered. To check the accuracy we should do error estimation to validate 
our results. If error is high, we have to go for adaptive mesh generation to reduce 
the error and to correct our solution to come closer to exact solution. Since we are 
in the nascent stage of this study, we are proceeding with our recovered smoothened 
stresses without doing any adaptive mesh generation using error estimation proce- 
dures. We are in a step prior to error estimation and in near future those modules 
will be added. 
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Chapter 5 


THE PROBLEM OF STIFFENER 
DESIGN 


Light weight vehicles(e.g. aeroplanes, light transport vehicles) have the outer 
structures made of aluminium alloys with suitable stiffeners to provide stiffness 
under bending and on shear loads. With this as a motivation, we consider the 
idealized sample stiffened plate configuration defined below. 


5.1 STIFFENED PLATE 


The stiffened laminated plate is taken as shown in Fig.5.1. The skin(top lamina) 
is taken to be made of aluminium alloy and the bottom stiffener is taken to be made 
of steel alloy. The goal being to remove material iteratively from the minimally 
stressed regions in stiffener such that the final configuration of the stiffened plate is 
an appropriately stressed minimum weight safe structure. The safe structure should 
satisfy the following conditions: 

1. The plate considered should not yield. 

2. The plate should not buckle. 
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Figure 5.1: Initial shape of stiffened plate 


3. The plate should not fail in vibrations. 

As an initial step in our study we are considering only the first condition. The 
procedure for optimizing stiffener topology is given below. 

5.2 PROCEDURE OF STIFFENER TOPOLOGY 
OPTIMIZAION 

A flowchart is given in figure 5.2 which shows the main stages in the stiffener 
topology optimization algorithm embedded in our program. User-defined parame- 
ters control the execution of the program and evolution of the structure. The details 
of the evolution are as given below. 

The initial structure is completely defined by the finite element mesh and this 
is read in together with applied loading and restraints. The finite element problem 
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Figure 5.2; Flowchart for stiffener design 
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4 divisions in x-axis and 4 divisions in y-axis 


Figure 5.3: Subelements in a master element 

is solved and smoothened stress field is recovered in postprocessing. The elements 
are then divided into subelements(see Fig.5.3) based on user specified divisions in x 
and y axis. The stress values at the vertices of the subelements are obtained using 




f \ 

(Jx 


^Xt 

<^y 

NN 

^yr 

Txy 

. = 2 jv; . 



1=1 


'Tyz 


'^yzi 

^XZ 


'^xzi 

V. > 

NN 



{<.} = -£Nl {cr.} (5.2) 

t=l 

where, NN is the number of nodes in an element, JV/ is the Lagrangian shape 
function associated with node in terms of normalised coordinates ^ and rj and 
{<Ti} is the generalized stress vector corresponding to node of an element. The 
element subdivision is carried in order to improve the accuracy of the identification 
of the low stressed material. 

A maximum cutoff value for stress is set for each lamina with respect to its 
corresponding yield stress. Von Mises stress is calculated at each subelement nodes 
using j 

(a) Mises = - (Tyf + 0-1 + al) + 34 + Zr^ + Zr ^'^ ' (5.3) 
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The acceptable design will be the one for which the Mises stress does not exceed 
the cut-off value at any point. Consider a lamina as shown in figure 5.4, For each 
node, the Mises stress is found at top, middle and bottom points of that lamina and 
the maximum value of those three points is stored for that node. If Mises stress 
values exceeds maximum cutoff value then stopping criterion is reached and the 
program is terminated. The finite element solution and mesh data of this initial 
configuration of stiffened plate is stored and evolutions are done using this data. 



Figure 5.4: A lamina with points defined 

For further evolution, the skin is untouched to maintain its external shape. Min- 
imum cutoflf value for stress is set with respect to yield stress for the stiffener such 
that the material in the regions having stress values less than this cutoff value is 
removed. Here the removal of material from an evolving design is typically done 
by rhanging the material properties, usually the Young’s modulus, to a near zero 
value, representing a softening of the material. The elements are said to be killed 
by decreasing the material property to near zero. Material is removed iteratively 
from the minimally stressed regions in the stiffener layer. In this way, in all regions 
of the design the stress will eventually be within a certain value of the maximum, 
producing a nearly fully stressed structure. To promote further change in the topol- 
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ogy of stiffener in each and every evolution, the minimum cutoff value for stress is 
increased incrementally based on a user supplied value. 

Stress contours are plotted for minimally and maximally stressed regions and re- 
gions are remeshed seperately and the procedure followed for plotting of stress con- 
tours and remeshing are given in the next section. Stress values at all the subelement 
nodes in the bottom laminate(stiffener) and the minimum cut off value are used to 
fit stress contours. The minimum stressed regions in stiffener are identified and the 
elements in that regions are ’killed’ by decreasing the material property to near zero 
and in the other regions material properties are unchanged. Then this evolved struc- 
ture is solved for and stopping criterion is checked. If stopping criterion is reached 
then further evolution is stopped else the finite element solutions and mesh data 
of our initial plate configuration is restored and minimum cutoff value for stress is 
increased. Then stress contours are plotted, regions are remeshed and the process is 
repeated. 


5.3 FITTING OF STRESS CONTOURS 

During the process of evolutionary optimization, at each step of the evolution 
stress contours, for a given stress level, are to be drawn in order to identify minimally 
and maximally stressed regions and then the different domains got after contouring 
are remeshed separately. Figure 5.5 shows a flowchart of the main stages in the 
process of fitting stress contours. The stages are listed below and described in detail 
below. 

Tracing the boundaries 
Polygons construction 
Splines construction 
Regions separation and Remeshing 
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Figure 5.5: Flowchart for fitting of stress contours 
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5.3.1 TRACING THE BOUNDARIES 


The initial step for fitting stress contours is to catch the boundary lines of 
minimally stressed regions in the stiffener. Steps followed for tracing the boundaries 
of minimally stressed regions are as follows 

1. Loop over the subelements of each and every element. 

2. Set the level of contouring for minimally stressed regions. 

3. Load the von Mises stress for the three corner nodes of that subelement. 

4. Checks are performed to find whether the subelement is having stress values, 
below the level of contouring or above the level of contouring or equal to the level 
of contouring or having stress values at some regions below the level of contouring 
and at some regions above the level of contouring. 

5. If the subelement has stress values below or above level of contouring then that 
subelements are checked for whether it is lying on the domain boundaries of our 
plate. If true, then that line (edge of that subelement) endpoints are noted. 

6. If the subelement has stress values equal to the level of contouring then the edge 
lines and endpoints of that element are noted. 

7. If the subelement has stress values which are above the level of contouring at 
some nodes and at other nodes it has stress values less than the level of contouring, 
then interpolation is done to find out line passing through the subelements which 
has stress values equal to that of level of contouring and the end points of that line 
is noted. 

8. After going through all elements we will get initial and final points of line segments 
of all the boundary contours of minimally stressed regions and also we will get the 
segmented lines of our plates domain boundary which are minimally and maximally 
stressed in separate files. 

5.3.2 POLYGON CONSTRUCTION 

Prom the continuity of stresses, there is a possibility of having only two types 
of contours (polygons) .They are closed and open contours. Closed contours will have 
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boundaries which will be totally inside the plate domain boundaries or it can start 
and end at the same point on a domain boundary of the plate. Open contours will 
have their end points on plates domain boundaries but at different locations. These 
two types of contours are separately identified. First we start with open contours, by 
identifying the end points on boundaries and start connecting the segmented lines 
in that element and then identifying its continuity in its neighbor element, like this 
search goes and connects each and every opencontours seperatly.Then the search 
jumps into connecting closed contours by starting with a line segment in a element 
in which the contour is passing and continuing in the same way till all the closed 
contours are separately connected. 

The boundary line segments are also connected to get boundary lines which are 
minimally and maximally stressed seperately.The following figure 5.6 shows how the 
boundaries are traced. 



Figure 5.6; Showing boundary plots ot various contours before spline construction 
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Figure 5.7: Showing boundary plots of various contours after spline construction 

5.3.3 SPLINE CONSTRUCTION 

Introduction 

The techniques for generating curves are of two types, the first one is the curve 
fitting technique. In this, the curve generation is constrained such that the curves 
pass through existing data points. This is particularly suited to shape description, 
where the basic shape is arrived at by experimental evaluation or mathematical 
calculations. Examples are aircraft wings, engine manifolds, mechanical and or 
structural parts. The alternative method uses free form of curves and surfaces 
which is suited for another class of shape design problems that depends on both 
aesthetic and functional requirements. Examples are ’skin’ of car bodies, aircraft 
fuselages, ship hulls, furnitures and glass ware. One form of generating free form of 

curves is B-splines. 

B-spline curves 

BYom a mathematical point ot view a curve generated by using the vertex of 
a dehning polygon is dependent on some interpolation or approximation scheme 
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to establish the relationship between the curve and the polygon. This scheme is 
provided by the choice of basis function. 

Here the B-spline basis (or polynomial approximation function) is generally non 
global. The nonglobal behaviour of B-spline curves is due to the fact that each 
vertex is associated with unique basis function. Thus each vertex affects the shape 
of a curve only over a range of parameter values where its associated basis function 
is nonzero. The B-spline basis (see [10]) also allows the order of the basis function 
and hence the degree of the resulting curves to be changed without changing the 
number of defining polygon vertices. 


The equation for B-spline curves can be expressed in a matrix form for both 
open and closed periodic B-spline curves. A open periodic B-spline curve is given 
by 


PAt’) = [r-]livi[Gl 

where here 

[T*] = .... 1] 0 < r < 1 

[G]’’ = [Bj B,+2 ... Bj+(ii-i)l 


l<i<n-l: + l, 0<r<l 


(5.4) 


j\r* 1 

•‘''i+l.m-l-l ~ {fc-l)! 


^ fc-1 ^ 


V 




m 


/ 


(5.5) 


\l-3 
0 < i,m < k — 1 


A closed periodic B-spline curve is given by 

Pj{t*) = [T*][iV*][G] 0 < y < n, 0 < t* < 1 (5-6) 


where here 

[T*] = .... 1] 0 < t* < 1 

[Gf = [B(j mod (n-H))+l ^(O'+l) rnod (n+l))+l - B(^^j+i+n-k) mod (n-f-l))-l-l] 


N; 




jfc-i 

i 


f 




l—m 


k 


0 < i,m < k — 1 


(5.7) 
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In the above expressions for open and closed B-spline curves Bj are the position 
vectors of the n+1 defining polygon vertices and k represents the order of the B- 
spline curve. Formally a B-spline curve is defined as a polynomial spline function of 
order k (degree k-1) since it satisfies the following two conditions: 

1. The function P(t*) is a polynomial of degree k-1 on each interval 0 < < 1 

2. p(t*) and its derivatives of order 1,2,3, ,k-2 are all continuous over the entire 

curve.Thus a fourth order is a piecewise cubic curve. 

Thus the output files from polygons construction which contains the vertex points 
of the polygons are fed into splines construction routine. Corresponding to open or 
closed polygon, open or closed splines are constructed. The degree of the curves 
generated are controlled by the user. Figure 5.7 shows the contours after spline 
construction. 

5.3.4 REGIONS SEPARATION AND REMESHING 

The open splines are connected to minimally stressed boundary lines to form 
closed minimally stressed regions,the same open splines are connected with maxi- 
mally stressed boundary lines to form closed maximally stressed regions.The closed 
contours are checked to find out which side of the closed contour is minimally 
stressed. Figure 5.8 shows how the regions are getting separated after spline con- 
struction. The main aim in regions separation is not to disturb the boundaries, if 
boundaries are distorted then we can’t proceed further for remeshing. 

Before proceeding to remeshing certain prerequisites are to be made to the output 
files which contains the seperated regions. 

1. The data contained in each file which represents the individual regions should be 
arranged in such a way that when a contour is drawn using the data from starting to 
ending, it should proceed in anticlockwise direction. Measures are taken to identify 
whether the given data points in a file is in anticlockwise direction or in clockwise 
direction. If the data points are in clockwise direction it is reversed and made 
into anticlockwise direction. The procedure followed for identifying whether data 
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Figure 5.8: Showing separated regions after spline construction 

points are arranged in anticlockwise direction or not is given after this prerequisites 
discussion. 

2. The files are to be arranged in such a way that the contours which are not 
enclosing closed contours are arranged first, then the contours which encloses closed 
contours are arranged next. 

3. The contours which encloses closed contours to be provided with informations 
that how many contours are inside that contours and also in which files are their 
informations are available. 

4. The data points in a file should be nearly equally spaced inorder to get elements 
of nearly equal sizes after remeshing. 

Figure 5.8 shows for the case taken how are the regions getting arranged and how 
is the direction of the vectors in the data files. Remeshing is done in each and every 
regions and it will always mesh to the left of a vector so that finally every regions are 
meshed inside. To assist this behaviour vectors are fed in anticlockwise directions. 
Figure 5.9 shows the remeshed regions. The procedure for finding whether a set of 
data points contained in a file are arranged in anticlockwise or in clockwise direction 
is as follows 
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Anticlockwise contour 


Clockwise contour 


Figure 5.10: Showing clockwise and anticlockwise contours 





1. The first point (a:i,yi) and second point (x 2 ,y 2 ) from the data file is taken, these 
points represents a line segment. A inward normal(P) which is patssing through 
midpoint of the linesegment is found using 

p ^ ^^1+^ _ 0.001(^2 - yi)) i + (p-~ + 0.001(12 - rri)) j (5.8) 

2. The inward normal(P) is shown in the figure for anticlockwise contour(5. 10(a)) 
and for clockwise contour(5. 10(b)). The point Xp and yp are given as 

_ o.001(y2 - yO) ; Vp = + 0.001(rr2 - rci)) (5.9) 

3. X-axis and Y-axis is taken at the point P{xp,yp) as shown in the figure. 

4. The angle included between the points 1 and 2 with respect to point P is found, 
similarly angles included between points 2 and 3 with respect to point P is found 
and gets added with initial one and this process continues till the angle included 
between last point and first point is added. 

5. The final included angle of the contour with respect to point F{xp,yp) turns out 
to be 360 degrees for anticlockwise contours and it is zero for clockwise contours. 
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Chapter 6 


NUMERICAL EXAMPLES 


6.1 PLATE GEOMETRY 


y- 



(Steel alloy) 


Figure 6.1: Stiffened plate geometry 

A laminated plate of the dimesions as shown in fig. (6.1) has been considered for 
analysis. 

The plate is loaded transversely on the upper surface with the sinusoidal load: 

q,{x, y) = -Qc sm(-Y-) sm(-^) (6.1) 
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The maximum cut-ofF value for stress is set as 80% of Yield stress of its corre- 
sponding material. The acceptable design will be the one for which the Mises stress 
does not exceed the maximum cut-off value at any point. 

The material properties for the top aluminium lamina are as follows: 

1. Yield stress(cry) = 406.81 10174A^/mm2 

2. Youngs modulus(E) = 72.3985709 x WN/mm? 

3. Rigidity modulus(G) = 27.58040796 x WN/mm?' 

4. Poissons ratio(//) = 0.33 

5. Density(p) = 2794.576 x l^~^Kg/mm? 

The material properties for the bottom steel lamina are as follows: 

1. Yield stress(<7j,) = 1213.53795iV/mm2 

2. Youngs modulus(E) = 199.9579 x lO^iV/mm^ 

3. Rigidity modulus(G) = 75.84612189 x lO^A^/mm^ 

4. Poissons ratio(i/) = 0.3 

5. Density(p) = 7830.34 x IQ-^Kglmm^ 

The boundary conditions used are simply supported and clamped edges. The de- 
grees of freedom fixed are shown in table 6.1 and 6.2. 


Simply Supported 

Clamped Edge 

o 

II 

II 

O 

il II 

lOo = 0 

Uq “ Ox ~ “ 0 

V0 = 6y=1py=0 


Table 6.1: Boundary Conditions on x edge ( y = constant) 


Simply Supported 

Clamped Edge 

1 

" II 
o 

II 

II 

o 

Wo = 0 

Uo = Ox = = Q 

Vo = Oy = i>y = 0 


Table 6.2: Boundary Conditions on y edge ( x = constant) 


44 









PROBLEM 1 

Four sides simply supported square laminated plate under sinusoidal transverse load 
on the top of the plate is analysed. The initial plate dimensions and load values are 
given below: 

X = b mm; Y = 5 mm; 
hi = 0.5 mm; h2 = 0.04 mm; 

Qc = 21N/mm^; m = 1 ; n = 1 
Shape functions of order 2 was used. 

Here the load taken is failure load for 0.5 mm thick aluminium plate without any 
stiffener and the problem is to find suitable stifiener design so that the plate does 
not yield and also stiffener is optimized to give minimum weight configuration. The 
thickness of steel plate is taken as 0.04 mm since at 0.03 mm the plate fails in the 
initial configuration itself. Here the material removal is done starting from minimum 
cutoff value of 10% of yieldstress of steel plate. At each iteration it is increased by 
1%. The plate is safe till the minimum cutoff value reached 33% of yieldstress of steel 
plate. It was noted that if the material removal is increased by increasing minimum 
cutoff value to 34% of yield stress of steel lamina, the Mises stress in the top lamina 
exceeded the maximu m cut-off value set for it. Since the problem is symmetric, 
quarter plate analysis was done and figures of quarter plate is given in figure 6.2. 
The figure 6.2(a) shows the top aluminium laminas stress profiles before material 
is removed from bottom steel lamina and figure 6.2(c) shows the stress profiles of 
top aluminium lamina after material was removed from the regions which are below 
or equal to 33% of yield stress of steel from the steel lamina. Figure 6.2(b) and 
6.2(d) shows the stress profiles of bottom steel lamina before and after material was 

removed. 

The area of stiffener before material was removed = 25mm 
The area of stiffener after material was removed = 16.8434017mm2 
Percentage of weight saved in stiffener = 32.626% 
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INITIAL CONF1GUR.AT10N 
(before material removal) 


Bottom Plate 

(Black regions are removed) 




(c) (d) 

WHITE REGIONS - Represents stress levels which are 65-80% of Yield stress of their corresponding material 
GREY REGIONS - Represents stress levels which are 33-65% of Yield stress of their corresponding material 
BLACk REGIONS - Represents stress levels which are 0-33% of Yield stress of their corresponding material 


Figure 6.2: Showing stress profiles for four sides simply supported square laminated 
plate (Quarter plate) 
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PROBLEM 2 

Four sides clamped square laminated plate under sinusoidal transverse load on the 
top of the plate is analysed. The initial plate dimensions and load values are given 
below: 

X = 5 mm; Y = h mm; 

hi = 0.5 mm; h2 = 0.0507 mm; 

Qc = 2ZN/mm^; m = 1 ; n = 1 

Shape functions of order 2 was used. 

Here the load taken is failure load for 0.5 mm thick aluminium plate without 
any stiffener and the problem is to find suitable stiffener design so that the plate 
does not yield and also stiffener is optimized to give miniTnnm weight configuratton. 
The thickness of steel plate is taken as 0.0507 mm since at 0.03 mm the plate fails 
in the initial configuration itself, and from 0.04 to 0.0506 mm the material removal 
is not significant. Here the material removal is done starting from minimum cutoff 
value of 5% of yieldstress of steel plate. At each iteration it is increased by 1%. 
The plate is safe till the minimum cutoff value reached 10% of yieldstress of steel 
plate. It was noted that if the material removal is increased by increasing minimum 
cutoff value to 11% of yield stress of steel lamina, the Mises stress in the top lamina 
exceeded the maximum cut-off value set for it. Since the problem is symmetric 
quarter plate analysis was done and figures of quarter plate is given in figure 6.3. 
The figure 6.3(a) shows the top aluminium laminas stress profiles before material 
is removed from bottom steel lamina and figure 6.3(c) shows the stress profiles of 
top aluminium lamina after material was removed from the regions which are below 
or equal to 10% of yield stress of steel from the steel lamina. Figure 6.3(b) and 
6.3(d) shows the stress profiles of bottom steel lamina before and after material was 
removed. 

The area of stiffener before material was removed = 25mm^ 

The area of stiffener after material was removed = 24.1316381mm^ 

Percentage of weight saved in stiffener = 3.473% 
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Top Plate 


INITIAI.CONFIGURVIION 
(before material remotal) 



m 


I 

f 

A 

I 





Bottom Plate 

(Black regions are removed) 



(a) (b) 

FINAL CONFIGURATION 
(after material was removed) 



(c) (d) 

WHITE REGIONS - Represents stress levels which are 65-80% of Yield stress of their corresponding material 
LIGHTGREY REGIONS - Represents stress levels which are 45-65% of Yield stress of their corresponding material 
GREY REGIONS - Represents stress levels which are 10-45% of Yield stress of their corresponding material 
BLACK REGIONS - Represents stress levels which are 0-10% of Yield stress of their corresponding material 


Figure 6.3: Showing stress profiles for four sides clamped squaie laminated 
plate{Quarter plate) 
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PROBLEM 3 

A square laminated plate with two sides simply supported( at X = 0 and X = 5) 
and other two sides free, under sinusoidal transverse load on the top of the plate is 
analysed. The initial plate dimensions and load values are given below; 

X = 5 mm] Y = 5 mm] 
hi = 0.5 mm; h2 = 0.01 mm; 

Qc = QN/mm^] m = 1 ; n = 1 
Shape functions of order 2 was used. 

Here the load taken is failure load for 0.5 mm thick aluminium plate without 
any stiffener and the problem is to find suitable stiffener design so that the plate 
does not yield and also stiffener is optimized to give min imum weight configuratton. 
The thickness of steel plate is taken as 0.01 mm . Here the material removal is 
done starting from minimum cutoff value of 10% of yieldstress of steel plate. At 
each iteration it is increased by 1%. The plate is safe till the minimum cutoff value 
reached 51% of yieldstress of steel plate. It was noted that if the material removal is 
increased by increasing minimum cutoff value to 52% of yield stress of steel lamina, 
the Mises stress in the top lamina exceeded the maximum cut-off value set for it. 
Since the problem is symmetric quarter plate analysis was done and figures of quar- 
ter plate is given in figure 6.4. The figure 6.4(a) shows the top aluminium laminas 
stress profiles before material is removed from bottom steel lamina and figure 6.4(c) 
shows the stress profiles of top aluminium lamina after material was removed from 
the regions which are below or equal to 51% of yield stress of steel from the steel 
lamina. Figure 6.4(b) and 6.4(d) shows the stress profiles of bottom steel lamina 
before and after material was removed. 

The area of stiffener before material was removed = 25mm^ 

The area of stiffener after material was removed = 9 . 20863876 Tn 7 n^ 

Percentage of weight saved in stiffener = 63.1654% 
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Top Plate 



INITIAL CONFIGUR.ATION 
(before material remo\ al) 




Bottom Plate 

(Black regions are removed) 


(a) (b) 

FINAL CONFIGURATION 
(after material was removed) 



WniTE REGIONS - Represents stress levels which are 65«80% of Yield stress of their corresponding material 
GREY REGIONS - Represents stress levels which are 51-65% of Yield stress of their corresponding material 
BLACk REGIONS - Represents stress levels which are 0-51% of Yield stress of their corresponding material 


Figure 6.4: Showing stress profiles for two sides simply supported and two sides free 
square laminated plate(Quarter plate) 


50 


Chapter 7 


CONCLUSIONS AND 
DISCUSSION 

7.1 CONCLUSIONS 

This study outlines an effective method to optimize stiffener topologies for plate 
structures. Prom the results obtained the following conclusions can be made 

1. Recovery algorithms gives better and smoothened stress field as compared to 
the one obtained directly from the finite element solution.This helps to improve the 
accuracy of results and also helps in plotting smooth stress contours, leading to 
smooth boundaries for the material regions. 

2. Capability to plot smooth contours of stress levels corresponding to a given 
threshold level using B-spline curves, removes the possibility of spurious corners be- 
ing formed after material is removed. 

3. Capability of remeshing removes the possibility of ‘wrong’ elements being re- 
moved. This remeshing algorithm also helps in changing the mesh size in each 
iteration according to users wish. 

4. The algorithm does not starts with a very fine mesh. The number of elements 
tends to increase during the evolution and hence it is computationally more efficient. 

5. The stiffener topologies depend on the type of boundary conditions and type of 


K.y.. -r-'-Tj’ 

V| , V I 1 

.; 5 , 
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loading conditions applied. Here it is shown that for change in boundary conditions 
not only the shape of the stiffener changes but also the thickness of stiffners also 
changes. Prom this, one should note that, to get a good design one should change 
both the thickness of the lamina and the stress levels. 

6. The algorithm proceeds in a semi-automatic way to help user to control the pro- 
cess according to the results obtained in the previous step. 


7.2 FUTURE WORK 

1. For a given mesh,element order p, laminate properties, transverse loads, our finite 
element program gives the solution corresponding to our plate model. However, in- 
order to make this code more effective the following modules need to be developed: 

(a) A-posteriori error estimation for both discretization and modelling error. 

(b) Adaptivity based on error estimation to get the optimal rate of convergence 
for the discretisation error, for any data of interest. 

(c) Curved element mapping to represent curved boundaries. 

(d) Refinement of mesh on the boundaries of contours, to minimize the discr- 
etization errors, instead of refining the mesh for the whole domain. 

2. Modules for applying different kinds of loading conditions can be added. 

3. Modules for buckling analysis, free vibration problems, etc can be added. 

4. The formulation can be generalized to curvilinear coordinates inorder to analyse 
shell structures. 

5. Code can be extended to include geometric non-linearity to handle large defor- 
mations. 

6. Instead of having single layer of stiffener, multilayered stiffeners can be optimized 
to get a 3-D profile for stiffeners. 

7. Modules for shape optimization can be added to obtain the optimized shape for 
the cutouts obtained from topology optimization. 

8. This method in principle is completely extendible to three dimensions without 
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any change except for the substitution of solid 3D elements for the 2D elements 
currently used. 
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APPENDIX A 
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